%% =========================== [Description] ============================
%% ---------------- implements the function -----------------
% Under multiple inputs
% Draw the mean variance curve of acceleration displacement         
% Draw a probability density curve for a given second             
% Draw response comparison curves for deterministic vs. random inputs
% ---------------- prefiles -----------------
% quasi0StiffSuspensionCertain.slx
% quasi0StiffSuspension.slx
% n1nDpolyRoots.m
% GetStatistics.m/GetStat.m
% ------------------ enter -------------------
% Cylinder air pressure
% Piston distance from back wall
%% ==================== [Select excitation signal and output] ========================
clear
disp('Pavement sine signal - input 1')
disp('Pavement random signal - input 2')
exc=input('Signal:');
disp('Output select acceleration - input 1')
disp('Output Select Displacement - input 2')
outputAmount=input('outputAmount:');
swiMat=[-1,1];
excSwitch=swiMat(exc);

%% ===================== [model run to find distribution] =============================
%% ----------------------------- [Determine values] ----------------------------------
d0Certain=32.5e-3; % piston distance from back wall m
mCertain=800; % mass on the spring kg
AcCertain=pi*0.1^2; % area of cylinder piston m^2
L0Certain=0.467; % cylinder node distance m
kCertain=3; %cylinder air pressure atm
m=mCertain;
m1=100; %Tire mass
Pa=101*10^3;
A_level=16; %-|
B_level=64; %-|-road
C_level=256; %-|
Level=B_level;
v=60; %Vehicle speed
C=2000; %damping   
order=3; 
D=1;
coeOfVar=[0,0,0,0.2];
%--------------------------------------------------------------------------
sim('quasi0StiffSuspensionCertain.slx')
Collocation=n1nDpolyRoots(order,D); 
groupNumber=length(Collocation(:,1));
opstable={accCertain,devCertain};
for groNum=0:groupNumber  
    if groNum==0
        [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                     coeOfVar,Collocation,groupNumber,1);
        sim('quasi0StiffSuspension.slx')
        outputData=zeros(length(acc.time),groupNumber);
    else
        [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                     coeOfVar,Collocation,groupNumber,groNum);
        sim('quasi0StiffSuspension.slx')
        switch outputAmount%Select the output amount
            case 1 
                output=acc.data;
                titleText='acceleration';
                yLabelText='Acceleration (m/s^2)';
            case 2
                output=dev.data;
                titleText='Displacement';
                yLabelText='displacement(m)';
            otherwise
                error('Select 1 or 2');
        end
        outputData(:,groNum)=output;
    end
   %-----------------------------------------------------------------------
    proPers=floor(groNum/groupNumber*1); %
    clc %---- progress bar area
    disp(['[','#'*ones(1,proPers),' '*ones(1,100-proPers) ...             %
        ,']',num2str(proPers),'%']) %
   % -----------------------------------------------------------------------
end
pointNumber=100;
timeNumber=10;
step=timeNumber/pointNumber;
time=0:step:timeNumber-step;
MCnumber=500;
PCmean=zeros(pointNumber,1);
PCvar=zeros(pointNumber,1);
MCdata=zeros(pointNumber,MCnumber^D);
for i=1:pointNumber    
% [PCmean(i),PCvar(i)]=GetStat ...
% (Collocation,outputData(100*i-99,:)',order); 
    [PCmean(i),PCvar(i),~,~,MCdata(i,:),~]=GetStatistics ...
    (Collocation,outputData(100*i-99,:)',order,MCnumber);
    %----------------------------------------------------------------------
    proPers=floor(i/pointNumber*99); %
    clc %---- progress bar area
    disp(['[','#'*ones(1,1),'#'*ones(1,proPers),' '*ones(1,99-proPers) ... %
        ,']',num2str(1+proPers),'%']) %
    % ----------------------------------------------------------------------
end
%% ====================== [Standard deviation and mean plot] ============================
figure(1)
errorbar(time,PCmean,sqrt(PCvar),'-o',...
    'color',[0,0,0],...
    'linewidth',0.8,...                        % Set line width     
    'MarkerSize',4,...                         % Set marker symbol size
    'MarkerEdgeColor','red',...                % Set marker edge size
    'MarkerFaceColor','red') %set the marker inner size
title([titleText,'Mean and standard deviation over time'])
xlabel('Time(s)')
ylabel(yLabelText)
disp('Program operation complete')
%% ==================== [Chaotic mean vs. model mean plot line] ====================
figure(2)
opDataCertain=opstable{outputAmount};
plot(time,PCmean,'-o','MarkerSize',4)
hold on
plot(opDataCertain.time,opDataCertain.data,'-.' ,'MarkerSize',4)
title([titleText,'Curve over time'])
ylabel('Output random variable mean')
xlabel('time(s)')
legend('random input','deterministic input')
%% ======================== [probability distribution for a given second] =========================
t=3.5;
tNumber=t/step+1;
[Xtpdf,Xt] = ksdensity(MCdata(tNumber,:));
figure(3)
plot(Xt,Xtpdf,'LineWidth',1)
title(['output',num2str(t),'probability density at s'])
xlabel('Output random variable')
ylabel('probability density \it{f}')
